Charge-4e superconductivity and chiral metal in 45°-twisted bilayer cuprates and related bilayers

The material realization of charge-4e/6e superconductivity (SC) is a big challenge. Here, we propose to realize charge-4e SC in maximally-twisted homobilayers, such as 45∘-twisted bilayer cuprates and 30∘-twisted bilayer graphene, referred to as twist-bilayer quasicrystals (TB-QC). When each monolayer hosts a pairing state with the largest pairing angular momentum, previous studies have found that the second-order interlayer Josephson coupling would drive chiral topological SC (TSC) in the TB-QC. Here we propose that, above the Tc of the chiral TSC, either charge-4e SC or chiral metal can arise as vestigial phases, depending on the ordering of the total- and relative-pairing-phase fields of the two layers. Based on a thorough symmetry analysis to get the low-energy effective Hamiltonian, we conduct a combined renormalization-group and Monte-Carlo study and obtain the phase diagram, which includes the charge-4e SC and chiral metal phases.

multi-component hexatic chiral superconductor leading to vestigial charge-6e SC was proposed in the context of kagome superconductors 22 .Presently, the material realization of the charge-4e/ 6e SC is still a big challenge.
Here in this work, we take advantage of the rapid development of the twistronics  and utilize it to design the intriguing charge-4e SC.
Here we shall study materials made through stacking two identical monolayers with the largest twist angle, which host Moireless quasicrystal (QC) structures [62][63][64] and are dubbed as the twist-bilayer QC (TB-QC) 65 , exampled by the recently synthesized 30 ∘ -twisted bilayer graphene [66][67][68][69][70] and 45°-twisted bilayer cuprates 71,72 .Prominently, the TB-QC hosts a doubly enlarged fold of rotation axis relative to its monolayer.Previous study 65,73 suggests that when each monolayer hosts a pairing state carrying the largest pairing angular momentum for the lattice, the second-order interlayer Josephson coupling (IJC) between the pairing ODPs from the two layers in the TB-QC makes them mix as 1: ± i, leading to time-reversal symmetry (TRS) breaking chiral topological SC (TSC).For example, as the monolayer cuprate carries the dwave pairing, the 45°-twisted bilayer cuprates will host the d + id chiral TSC 65,[74][75][76][77] .It's interesting to investigate possible vestigial secondary orders above the T c of these chiral TSC phases, driven by the secondorder IJC between the pairing ODPs from the two layers.
In this paper, we study the secondary orders in the superconducting TB-QC.Its unique symmetry leads to a simplified lowenergy effective Hamiltonian including decoupled total-and relativephase fields between the bilayer.Significantly, the second-order IJC allows the relative phase to fluctuate between its two saddle points to restore the TRS.Consequently, while the unilateral order of the relative-phase field leads to the TRS-breaking chiral-metal phase, the unilateral quasi-order of the total-phase field leads to the charge-4e SC phase, in which two Cooper pairs from different layers pair to form quartets.These two vestigial phases occupy different regimes in the phase diagram obtained by our combined renormalization group (RG) and Monte-Carlo (MC) studies, which are unambiguously identified by various temperature-dependent quantities including the specific heat, the secondary ODPs, and their susceptibilities, as well as the spatialdependent correlation functions.

Model and symmetry
Taking two D n -symmetric monolayers, let's stack them by the twist angle π/n to form a TB-QC, as shown in Fig. 1 for n = 4 (e.g., the cuprates) and n = 6 (e.g., the graphene).Obviously, the point group is D nd , isomorphic to D 2n .There is an additional symmetry generator in the TB-QC which is absent in its monolayer, i.e., the C 1 2n rotation accompanied by a succeeding layer exchange, renamed as C1 2n here.Suppose that driven by some pairing mechanism, the monolayer μ = t/b (top/bottom) can host a pairing state with pairing angular momentum L = n/2.While the cuprate monolayer hosting the d-wave SC synthesized recently 78 provides a good example for n = 4, some members in the graphene family which were predicted to host the f-wave SC 73,79,80 set an example for n = 6.The pairing gap function in the μ layer is Δ ðμÞ ðkÞ = ψ μ Γ ðμÞ ðkÞ: ð1Þ Here Γ (μ) (k) is the normalized real form factor, and ψ μ is the "complex pairing amplitude".Prominently, the Γ (μ) (k) for L = n/2 changes sign with every C 1 n rotation.As shown in Fig. 1, we choose a gauge so that Γ ðbÞ ðkÞ = Pπ n Γ ðtÞ ðkÞ, P2π n Γ ðμÞ ðkÞ = À Γ ðμÞ ðkÞ: ð2Þ Here Pϕ indicates the rotation by the angle ϕ.As the interlayer coupling in the TB-QC is weak 62,63,65 , we can only consider the dominant intralayer pairing, but the two intralayer pairing ODPs can couple through the IJC [73][74][75][76][77]80 . We hall investigate the ground state and the vestigial secondary orders induced by this IJC.
Firstly, let's make a saddle-point analysis for the Ginzburg-Landau (G-L) free energy F as functional of ψ t/b .For the saddle-point solution, the ψ t/b are spatially uniform constant numbers.F is decomposed as, where F 0 ðjψ μ j 2 Þ are the monolayers terms and F J is the IJC.The TRSallowed first-order IJC takes the form, F Under C1 2n , the gap function on the μ layer changes from Δ ðμÞ ðkÞ = ψ μ Γ μ ð Þ ðkÞ to ΔðμÞ ðkÞ = ψ μ Pπ n Γ μ ð Þ ðkÞ which, under Eq.( 2), can be rewritten as ψμ Γ μ ð Þ ðkÞ with The invariance of F under C1 2n requires α = 0. Thus, the following second-order IJC should be considered, Equation ( 6) is minimized at ψ b = ±iψ t for A 0 > 0 or ψ b = ±ψ t for A 0 < 0. Previous microscopic calculations favor the former for the 45°twisted bilayer cuprates 65,74 and 30°-twisted bilayer of the graphene family 73,80 , leading to d + id or f + if chiral TSCs ground state.
Secondly, let us provide the low-energy effective Hamiltonian for the pairing-phase fluctuations.In this study, we fix Γ (μ) and set ψ μ !ψ μ r ð Þ as a slowly varying "envelope" function to describe the spatial fluctuation of the complex pairing amplitude.Focusing on the phase fluctuation, ψ t/b are written as ψ t=b = ψ 0 e iθ t=b r ð Þ where ψ 0 > 0 is a constant.The θ t=b r ð Þ are further written as Here θ + r ð Þ and θ À r ð Þ denote the total and relative pairing phases.The low-energy effective Hamiltonian reads with ∂ ± ≡ ∂ x ± i∂ y .Up to the lowest-order expansion, the H 0 takes the following explicit form in the k-space, 2n , the gap function on the μ layer changes from 2), can be rewritten as Consequently, we have The invariance of Eq. ( 9) under ( 11) only allows for nonzero ρ and κ, leading to the real-space Hamiltonian Equation ( 12) shows two important features.Firstly, the θ + and θ − fields are dynamically decoupled, with each hosting different stiffness parameter ρ or κ derived by the G-L expansion in the Sec.I of Supplementary Information (SI).Secondly, the second-order IJC allows θ t − θ b = 2θ − to fluctuate between its two saddle points, i.e., ±π/2, to restore the TRS.Note that although the term cosð4θ À Þ in Eq. ( 12) leads to four different values of θ − : ±π/4 and ± 3π/4 for the ground state, θ − = π/4 (−π/4) leads to gauge equivalent state with θ − = −3π/4 (3π/4).So the system only possesses two-fold Ising anisotropy.Here the unilateral quasi-ordering of the θ + field leads to the ODP Δ (t) (k) ⋅ Δ (b) (−k) characterizing the charge-4e SC in which two Cooper pairs from different layers pair.The unilateral ordering of the θ − field leads to the ODP Δ (t)* (k) ⋅ Δ (b) (k) characterizing the TRS breaking chiral metal 81,82 .Note that while θ + and θ − each can host either integer or half-integer vortices, Eq. ( 7) requires that they can only simultaneously host integer or half-integer vortices to ensure the single-valuedness of ψ t/b 10,18 . This sets the "kinematic constraint" in the low-energy "classical Hilbert space" for allowed vortices of the two fields.
RG study: To perform the RG study, we start with the following effective action at the temperature T, Here g 4 > 0 is proportional to A 0 .This action can be mapped to a two-component Sine-Gordon model, The dual bosonic fields θ + and θÀ describe the vortices of the fields θ + and θ − .g 2,0 , g 0,2 , and g 1,1 are coupling parameters proportional to the fugacities of different types of vortices (g 2,0 /g 0,2 : integer vortices; g 1,1 : half vortices).
The phase diagram obtained by the one-loop RG analysis provided in "Methods" is shown in Fig. 2a.Variation of the initial coupling parameters does not change the topology of the phase diagram, which always includes the chiral TSC, charge-4e SC, chiral metal, and normal metal phases, see the Section IV of SI.At low enough T, the vortex fugacities g 2,0 , g 0,2 , and g 1,1 are all irrelevant while the IJC parameter g 4 is relevant, suggesting that both the θ ± fields are locked, leading to the TRS breaking chiral SC.With the enhancement of T, in the low κ/ρ regime, g 0,2 first gets relevant (and suppresses g 4 ) suggesting that the θ − vortices proliferate to restore the TRS, to form the charge-4e SC.In the high κ/ρ regime instead, the g 2,0 first gets relevant suggesting the θ + vortices proliferate to kill the SC, to form the chiral metal.In both regimes, at high enough T, g 2,0 , and g 0,2 are both relevant, forming the normal metal phase.In the regime κ ≈ ρ, with the enhancement of T, the system transits into a phase wherein the coupling g 1,1 is relevant and the half vortices involving both fields proliferate to kill both (quasi) orders, suggesting that the system directly transit to the normal state.
In the charge-4e SC, the Josephson-coupling phase, i.e., θ − , is disordered.However, this phase should not be understood as a layerdecoupled charge-2e SC from each layer, as in this phase the pairing phase of each layer is also disordered.To remind you, the charge-4e SC proposed here only lives in the intermediate temperature above the T c of the pairing state, wherein each layer is no longer superconducting.In the chiral metal phase, the time-reversal symmetry breaking can be verified by the polar Kerr effect.Furthermore, there can be spontaneously generated inner magnetic field in the material, which can be detected by the muon spin resonance experiment.

MC study
To perform the MC study, we discretize the Hamiltonian (12) on the square lattice to obtain Here 〈ij〉 represents nearest-neighbor bonding, and the positive coefficients α, λ, and γ satisfy, Note that although different α, λ, and γ satisfying Eq. ( 16) leads to the same continuous Hamiltonian (12) in the continuum limit, it is required that all of them should be positive so as to reproduce the correct low-energy "classical Hilbert space" for allowed vortices.The reason is as follows.Here the α > 0 and λ > 0 terms energetically allow for integer or half-integer θ + and θ − vortices, while the γ > 0 term energetically only allows for integer θ t or θ b vortices and hence imposes the "kinematic constraint" between the θ + and θ − vortices.Note that although the γ term does not naturally emerge from Eq. ( 12), the singlevaluedness of the ψ t/b field dictates it.This term is crucial to yield the correct topology of the phase diagram.As shown in the Sec.VI of SI, if we turn off the γ term, θ + and θ − are decoupled, leading to a topologically wrong phase diagram.A comparison between the correct phase diagram and the wrong one shows that the kinematic correlation makes the vestigial phase regimes largely shrink.For thermodynamic limit, even an infinitesimal γ can energetically guarantee the "kinematic constraint".Here in the discrete lattice, we set γ = 1 4 ρκ=ðρ + κÞ, A = 0:025ρ, and their other values lead to similar results, see Section VI of SI.
The MC phase diagram shown in Fig. 2b is qualitatively consistent with the RG one shown in Fig. 2a.Various T-dependent quantities are shown in Fig. 3 for κ/ρ = 0.3, 1, 2.2 marked in Fig. 2b, with the formulas adopted in the MC calculations provided in "Methods".For κ/ρ = 0.3, the specific heat C v is shown in Fig. 3a, where the high-T broad hump characterizes the Kosterlitz-Thouless (K-T) phase transition between the normal state and the charge-4e SC and the low-T sharp peak characterizes the Ising phase transition between the charge-4e SC and the chiral SC.For this κ/ρ, Fig. 3b shows the phase stiffness S characterizing the SC and the Ising ODP I characterizing the relative-phase order 16 , which emerge at the critical temperatures corresponding to the broad hump and sharp peak in Fig. 3a, respectively.Furthermore, the total-(χ + ) and relative-(χ − ) phase susceptibilities 16 shown in Fig. 3c diverge at the same critical temperatures.For κ/ρ = 1, the specific heat shown in Fig. 3d exhibits only one peak, suggesting a direct phase transition from the normal state to the chiral SC.Such a result is also reflected in Fig. 3e, f which shows that the total-and relative-phase (quasi) orders emerge at the same temperature.For κ/ρ = 2.2, the corresponding results shown in Fig. 3g-i reveal that following the decrease of T, the system will successively experience the normal state, the chiral metal, and the chiral TSC phases.The results presented in Fig. 3 are well consistent with the phase diagram shown in Fig. 2b.
The total-(+) and relative-(−) phase correlation functions η ± are shown in Fig. 4. See their formulas in "Methods".Figure 4a, b shows that for the representative point A marked in Fig. 2b, while η + (Δr) power-law decays with Δr suggesting quasi-long-range order of the total phase, η − (Δr) decays exponentially with Δr, suggesting disorder of the relative phase.Obviously, these electron correlations are consistent with the charge-4e SC phase.Figure 4c, d shows that for the point D, while η + (Δr) decays exponentially with Δr suggesting disorder of the total phase, η − (Δr) saturates to a constant number for large enough Δr suggesting long-range order of the relative phase, consistent with the chiral-metal phase.For comparison, the η ± for points B and C provided in Section V of SI is also consistent with the normalmetal and chiral-SC phases.

Discussions
In comparison with previous proposals for the charge-4e/6e SC based on melting of the PDW 10,11,14 or the nematic pairing 17,18 , our proposal is based on a more definite and easily realized start point: here we only need to start from non-topological d-wave SC (or f-wave SC) in any fourfold (or sixfold) symmetric monolayers.Particularly, we have provided concrete synthesized materials to realize our proposal, i.e., the 45 o -twisted bilayer cuprates and the 30 o -twisted bilayer of some graphene family.Furthermore, superior to the previous bilayer approach, here a Cooper pair from the top layer pairs with a Cooper pair from the bottom layer to form the charge-4e SC between the layers.Consequently, the half flux quantization can be experimentally detected as a hallmark of the charge-4e SC in our proposal.
The TB-QC provides a better platform to realize the vestigial phases than conventional chiral superconductors such as the p + ip or d + id ones on the square or honeycomb lattices.The latter also hosts two degenerate pairing ODPs and hence can accommodate both total and relative-phase fluctuations of the two ODPs.However, the rotational symmetry of the monolayer system is not as high as that of the TB-QC studied here.Consequently, for chiral TSC in monolayer systems, there can be many nonzero coefficients in Eq. (9).Particularly, the two-phase fields are generally dynamically coupled as the symmetries in these systems allow for extra terms such as ∇ ± θ + ⋅ ∇ ± θ − in the Hamiltonian density in Eq. (12).See more details in Section II of SI.As shown in Fig. 2 and Fig. S5a, the kinematic correlation between θ + and θ − has already made the vestigial phase regimes largely shrink, their extra dynamic coupling might make them further shrink or even vanish.
In conclusion, we have predicted the realization of the charge-4e SC or the chiral metal in the TB-QC, emerging as the unilateral (quasi) ordering of the total-or relative-pairing phase of the two layers, above the chiral-TSC ground state.The TB-QC provides a better platform to realize these vestigial phases than previous proposals as here we can start from a more definite and easily realized start point.

Methods
The RG approach Here we provide some technique details for the RG study.With standard RG analysis, the flow equations at the one-loop level are given by: Here b represents the renormalization scale, g 2,0 , g 0,2 , and g 1,1 represent the coupling strength of different types of topological defects, ρ 0 = ρ=T and κ 0 = κ=T represent two kinds of stiffness parameters.
The fixed points of N general RG flow equation dg d' = RðgÞ is obtained by R(g * ) = 0.In Table 1, we present four fixed points of the RG flow equation ( 17) and the corresponding phases in our calculation.We find the renormalized values of the stiffness parameters ρ 0 and κ 0 are consistent with the phase revealed by the RG flow result of the g-couplings.Specifically, the ρ 0 flows to a finite positive value if the U(1)-gauge symmetry is (quasi) broken, otherwise it flows to zero; the κ 0 flows to infinity if the time-reverse symmetry is broken, otherwise it flows to zero.Furthermore, the stability analysis of the fixed points can be provided following the standard process 83 .We only outline here and see more details in the Section III of SI.The β function of the coupling constant which is very close to the fixed point g * can be replaced by a linear mapping: where we have used R(g * ) = 0, and M αβ = ∂R α ∂g β j g = g * .To get the stability properties of the flow, we have to diagonalize the matrix M N×N .The eigenvalues denoted by λ α , α = 1, 2, . . ., N. If the real parts of all the eigenvalues are negative or, at worst, zero, i.e., the scaling fields are all irrelevant or marginal.There are stable fixed points corresponding to the "stable phases".Complementary to the stable fixed points, if all the eigenvalues are positive and the scaling fields are all relevant, there are unstable fixed points.Additionally, there is a generic class of fixed points with both relevant and irrelevant scaling fields.These points are associated with the boundary of the phase transition.

The Monte-Carlo approach
Here we provide some formulas for the MC calculations.The phase stiffness characterizes the quasi-long-range order of the total phase and hence the SC is 16

S
Table 1 | Fixed points of the coupling parameters under RG, and the corresponding phases with where N is the site number, and β = 1/k B T.
The Ising order parameter characterizing the relative-phase ordering breaking the time-reversal symmetry is, The total-(+) and relative-(−) phase susceptibilities for temperatures above the T c of the corresponding orders are defined by The total-(+) and relative-(−) phase correlation functions are defined as X r e i½θ t ðrÞ ± θ b ðrÞÀθ t ðr + ΔrÞ∓θ b ðr + ΔrÞ

Fig. 1 |
Fig. 1 | Schematic illustration of a TB-QC formed by two D n -symmetric monolayers, with each monolayer carrying SC with pairing angular momentum L = n 2 .As examples, n = 4 (cuprates) and n = 6 (graphene) are shown in (a) and (b), respectively.

Fig. 2 |
Fig. 2 | The phase diagram.Phase diagram provided by a the RG study and b the MC study.The initial values of the coupling parameters in a are g 2,0 = g 0,2 = 0.1, g 1,1 = g 4 = 0.01 in Eq. (14) and in (b) are A = 0.025ρ and γ = 1 4 ρκ=ðρ + κÞ in Eq. (15).The three dotted lines and four dots (A-D) in (b) are used for subsequent interpretation of the phase diagram.

Fig. 4 |
Fig. 4 | The correlation functions.The correlation function η ± for a and b for A(κ = 0.2, T = 0.2), and for c and d for D(κ = 2.2, T = 0.5) marked in Fig. 2b.Insets of a the log-log plot, and b and c only the y-axes are logarithmic.